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ABSTRACT 

We investigate the chaotic behavior of the S-star cluster in the Galactic center using precise N-body calculations¢ free from 
round-off or discretization errors. Our findings reveal that chaos among the Galactic center S-stars arises from close encounters, 
particularly among pairs and near the massive central body. These encounters induce perturbations, causing sudden changes in the 
orbital energies of the interacting stars. Consequently, neighboring solutions experience roughly exponential growth inseparation. 
We propose a theory of "punctuated chaos" that describes the S-star cluster’s chaotic behavior. This phenomenon results from 
nearly linear growth in the separation between neighboring orbits after repeated finite perturbations. Each participating star’s orbit 
experiences discrete, abrupt changes in energy due to the perturbations. The cumulative effect of these €vents is further amplified 
by the steady drift in orbital phase. In the Galactic center, perturbations originate from coincidental encounters occurring within 
a distance of < 100 au between at least two stars (in some cases, three stars). Our model satisfactorily explains the observed 
exponential growth in the 27 S-star cluster. We determine that the S-star system has a Lyapunov time scale of approximately 
462 + 74 years. For the coming millennium, chaos in the S-star cluster will be driven.mainly by a few of the closest orbiting 
stars: S2, S5, S6, S8, S9, $14, $18, S31, $21, $24, $27, S29, and S38. 
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1 INTRODUCTION onant relaxation in the galactic center (Sridhar & Touma 2016), or 
even with violent relaxation (Kandrup et al. 2003). Previous studies, 
however, found that repreated weak scatterings among minor bod- 
ies also form a major driver for changes in the orbital paramaters, 
for example in studies on the orbit of mercury (Laskar & Gastineau 
2009), but also in Halley’s comet (Boekholt et al. 2016) to be chaotic 
on a time scale of less than 3000 years, as it interacts with Venus 
and Jupiter. In numerical studies of planet-planet scattering experi- 
ments in hypothetical multi-planet systems chaotic behavior is also 
recognized to be driven by encounters, rather than resonant overlap 
(Jurić & Tremaine 2008; Chatterjee et al. 2008). This event-driven 
chaotic behavior does not comply with Chirikov’s resonant overlap 
paradigm. 


Newton’s laws of motion lead to chaos. This chaotic behavior is 
often quantified by measuring the growth in time of the,separation, 
ô, between two neighbouring solutions, i.e. by solvingthe equations 
of motion of the multi-body system twice, once withrandonce without 
a small change in the initial conditions. If the evolution of 6 is roughly 
exponential, but 6 is still small by the end of,the calculations, we call 
the evolution chaotic. 

Chaotic behavior can be quantified.using the Lyapunov time scale, 
which (in a more rigorous treatment) is, the reciprocal of the maxi- 
mum positive Lyapunov exponent. The Lyapunov time scale has been 
measured for only a limited number of multi-body systems, because 
these measurements aresexpensive in terms of computer time. It is 
not even clear that a(specific system, within the uncertainty of its 
initial conditions/leads to_atinique Lyapunov timescale, because the 
available parameter space may have an irregular structure with stable 


This finding hints towards an event-driven process, where reso- 
nances are not required, but where chaos is driven by close encoun- 


and chaotic Tegions.¢Hayes 2008). 

In the’Solar System, chaos is mainly driven by resonant overlap 
(Chirikovel9797 Wisdom 1980; Tamayo et al. 2021; Rath et al. 2022; 
Mogavero & Laskar 2022), which has some resemblance with res- 
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ters. We informally refer to such behaviour as “punctuated chaos”. 
A similar chaotic behavior was found in stellar clusters (Goodman 
et al. 1993), where mutual interactions between stars seemed to drive 
chaos in the system. A system consisting of many eccentric and in- 
clined orbits behaves differently from a flat system with non-crossing 
orbits and a dominant central body. The S-star cluster in the Galac- 
tic center (Ghez et al. 2008; Gillessen et al. 2009) may be a good 
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example where resonant overlap does not form the major driver for 
chaos. 

Understanding chaos in self-gravitating systems is important for 
understanding a wide variety of astronomical phenomena, including 
the sources of gravitational waves, extreme-mass-ratio inspirals, and 
the probability that a minor body hits the planet Earth. Our picture 
of punctuated chaos is also important for other Hamiltonian systems, 
such as the multi-body pendulum (Hoover & Griswold Hoover 2017), 
and it may be responsible for a slow-down of chaotic behavior in 
critical fluctuations (Das & Green 2019). 

In our analysis of the S-stars, we will neglect low-mass stars, 
unobserved compact objects, and other interstellar material (planets, 
asteroids, dust, gas). The low-mass components are expected to be 
copious and important, but we limit ourselves to study chaos in an 
idealized system of black hole with S-stars. Chaos in the here studied 
idealised system then probably only provides an upper limit to the 
actual Lyapunov time scale. 


2 CHAOS AS AN EVENT-DRIVEN PROCESS 
2.1 The Lyapunov timescale 


Consider a Keplerian two-body system, e.g. a star with a planet or 
supermassive black hole with a star. The difference between two 
neighbouring solutions, started with slightly different initial condi- 
tions, grows approximately linearly with time. The orbital phases 
of the two solutions gradually diverge, until their separation 6 has 
grown to the size of the orbits. 

In a three-body system consisting of a binary and a third body, the 
binary orbits are repeatedly perturbed by the third body. The effect 
of these events on two neighbouring binary orbits depends on their 
separation 6, and can cause the two solutions to diverge more quickly 
than in a two-body system. A sequence of perturbations, each with 
a subsequent period of linear divergence until the next perturbation, 
can then drive exponential divergence between two initially almost 
identical systems (Sec.2.2). The degree of chaos, measured in terms 
of the Lyapunov time scale, can then be derived by the frequency and 
strength of perturbations in the system (Sec.4. 1). 

In direct N-body calculations we can quantify chaos by measuring 
the Lyapunov time scale. Instead of a measure over the system’s 
entire lifetime, which would be the setting for a rigorousAreatment/a 
measure over a finite time interval is more practical in our’case. We 
therefore define the growth factor G5(t) following from an initial 
separation 6(0) after some time t. The evolution of the separation is 
then described as an exponential function of time 5(t) = 6(O)e*, 
where J is the Lyapunov exponent. The growth factor can then be 
written as G 5 = 6(t)/d(0) = e% , and A = log(@s)/t, or in terms of 
the Lyapunov timescale, t4 = 1/4 = A/log(G>5). 


2.2 Punctuated chaos 


Based on the event-driven process described above, we derive ap- 
proximate expressions for the Lyapunov time scale. We start again 
by considering a particle in a Kepler orbit around a much more mas- 
sive body. With‘an initial semi-major axis ag, and total mass m the 


initial orbital frequency wo = ,./Gm/ ap. Let a neighboring solution 
be separated by’an infinitesimal displacement ôxọ at time t = 0. This 


displacement leads to a small difference in the semi-major axis of the 
same order, i.e. dag ~ Oxo. The resulting difference in frequency is 


ôwo ~ 5x9fGm/ap. (1) 


The separation of the two motions along the orbit grows with time 
t > 0 according to 


6x(t) ~ 6x9 + bwWoaot = 6xg(1 + wot). (2) 


The growth of the initial displacement is linear with time from tọ to 
t, but such that 6a remains constant as the growth is along the orbit, 
i.e. the growth is in orbital phase rather than in energy (or angular 
momentum). 

Now suppose that an instantaneous perturbation acts on the motion 
at time rt, , causing the velocity of the Kepler motion to receive a slight 
kick. 

There is a resulting change in energy, and therefore a change in 
semi-major axis Aa]. But since the perturbation will depend on the 
position of the body, the difference in semi-major axis between the 
two motions, da, will also change, by an amount which we call 
ôa. We suppose for the sake of argument that da; ~ 6x,,A.e. the 
difference in position at the time of the event. (In Sec.4. lawe give a 
specific example in which it is easy to see that da; œ ôf.) The new 
difference in semi-major axis implies a difference in orbital frequency 
Ow, ~ 6a,w ,/a, at time tı. Thus for t > tı the displacement varies 
as 


ôx ~ 6x, + 1(t-t)da, 
~ 6x11 +a) (at). (3) 


If a second perturbation occurs at times) > t], we can see from 
eqs.2 and 3 that the displacement is 


0x2 ~ dxo FF (1 + w(t - t1)). (4) 


If these perturbations\recur aðfoughly comparable intervals At, and 
if w does not change by aarge factor, it can be seen that the dis- 
placement at some largetime ¢ will be 


x(t) ~ sxo + wA "A. (5) 


Thus‘the linear growth of eq. 2 transforms into exponential growth. 
The corresponding Lyapunov exponent is O(w) if wAt < 1; it is of 
order the reciprocal of the crossing time. The case wAt Z 1 is also of 
interest and leads to a smaller estimate of order In(wAr) /At. 
Though we here considered perturbed Keplerian motion with one 
dominant body, which is relevant to our discussion of the S-stars 
in Sections 3 and 4, the conclusions are valid more widely. As an 
example, we consider resonant three-body scattering events, which 
can be viewed as a prolonged sequence of perturbations of the Kepler 
motion. So long as the three particles remain democratic (i.e., at 
comparable distances) and are of comparable mass, the perturbations 
in any of the three two-body motions will be of order unity and 
occur at intervals of order the crossing time, ter. Therefore, as in 
the above discussion, the Lyapunov exponent will be of order 1/ter. 
Numerical examples of Dejonghe & Hut (1986) indeed show that 
the separation of neighboring solutions grows roughly exponentially 
until dissolution of the resonance (in the sense of a democratic three- 
body behaviour). On the other hand, if the evolution of a triple system 
is dominated by protracted excursions of the third body, of order T > 
tcr, then the estimate will decrease to one of order 1/T (following the 
result for the case wAr Z 1, and neglecting a logarithm). Usually, the 
evolution is a mix of prolonged excursions interspersed with periods 
of frequent interplay (Szebehely 1971), and the Lyapunov exponent 
A will be 1/T $ A < 1/ter, where T is the duration of the longest 
excursion. The lifetime of the Pythagorean problem, for example, 
is about 16fcr (Aarseth 2003, their p.238), and the growth of the 
separation of neighboring solutions in this time is about 8.5 dex 
(Dejonghe & Hut 1986). Thus T $ 16tcer and A = 1.2/ter, and the 
finite-time Lyapunov exponent lies in the range 1/T $ A £$ 1/ter. 
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The result of the model (that the Lyapunov exponent J is of order 
1/tcr for comparable masses without lengthy excursions) is consistent 
with the results in Goodman et al. (1993), who considered the gen- 
eral N-body problem. This is a rather independent confirmation, as 
their model was also based on assuming that the separation between 
neighboring solutions grows as a result of two-body encounters (see 
Tab. 1, but otherwise did not follow the same approach as ours. 


3 APPLICATION OF PUNCTUATED CHAOS 


We now study punctuated chaos in an actual astrophysical system, 
namely the S-star cluster in the Galactic center. All stars affect each 
other and we therefore stretch the assumptions of extreme mass ra- 
tio, the approximation of only two orbiting particles, and the three- 
dimensional nature of the problem. 

Our analysis starts by acquiring converged solutions for the dy- 
namical evolution of the S-stars. Acquiring a converged solution is 
important because round-off and discretization errors grow exponen- 
tially, rendering a non-converged solution inappropriate for studying 
chaos. We acquire converged solutions to the N-body problem by 
integrating Newton’s equations of motion. The dynamics near super- 
massive black holes is best described with general relativity rather 
than by Newtonian dynamics. Still we perform our simulations using 
a Newtonian approximation and ignore the theory of general rela- 
tivity. We consider this appropriate because the highest speed in our 
simulations (the star S2 at pericenter) does not exceed about ~ 0.017 
times the speed of light. At these relatively low speeds, the system 
behaves like a Newtonian system, at least within our brief simula- 
tion time of 104 years (Portegies Zwart et al. 2022). In addition, our 
model for punctuated chaos was derived for Newtonian systems. In 
the presence of the black hole, in future it will be worth exploring 
the effect of general relativity on the chaotic behavior of the S-star 
cluster. 

The results presented here are acquired using converged solutions. 
These result from a procedure in which the length of the mantissa, 
controlling precision, and the accuracy of the integrator are improved 
at each iteration step, as is explained in Portegies Zwart & Boekholt 
(2018). Such converged solutions can be achieved using Brutus: In 
Brutus, we control round-off by extending the numerical preCision, 
and accuracy by reducing the tolerance in the Bulirsch-Stoer integra- 
tor (Bulirsch & Stoer 1964). By repeating the same caléulation with 
higher precision and better accuracy, we eventually reach a solution 
for which the results remain identical to a pre-determined,number of 
decimal places; we call this the converged solution) 


3.1 Measuring chaos in the S-star cluster 


We consider a realization of the S-star cluster, consisting of 27 early- 
type stars that orbit the supermassive black hole in the Galactic 
center. We adopted the orbital parameters of the 27 S-stars reported 
in Gillessen et al. (2009)(their.table 7). The numbering of these stars 
is not simply S1 to $27, because we only use those stars for which 
an orbit is determined. 

The initial, conditions ‘are generated using the AMUSE (Portegies 
Zwart et“ “al.2013; Pelupessy et al. 2013) routine 
generate_Sstar_cluster.py from Portegies Zwart & McMillan 
(2018). This routine adopts the cited orbital elements for the S-stars, 
and solves’’Kepler’s equation to acquire Cartesian coordinates 
and velocities of the S-stars. We assume each S-star to have a 
mass of 20 Mo, and adopt the central massive black hole mass of 
4.154 - 10° Mo (Gravity Collaboration et al. 2019). We integrate 
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Figure 1. Projection of the S-stars’ orbits in the y-x plane. These orbits 
are integrated for10.000 years from January 1st 2001 to 12 001 yr. Even the 
widest orbits (star$83 has an orbital period of ~ 1700 yr) are overplotted 
several times) The orbits look Keplerian, but when we zoom in on any orbital 
segment, the fine structure of the chaotic orbital evolution becomes visible 
(seefigure’2). 


this system for 10000 yr until the solution converges using Brutus 
(Boekholt & Portegies Zwart 2015). 

In figure 1 we present the converged solution of the projected 
orbits of the S-stars integrated over 10000 years. This solution was 
obtained using a word-length Ly = 128 bits, and a tolerance in 
the Bulirsch-Stoer integrator of € = 10-24 leading to a total final 
phase-space error less than 1/107, and a relative energy error in the 
integration of dE/E ~ 107!5. Each calculation took about 20 hours 
on the single core of a Xeon W-11855M-processor. For Ly = 116 bit 
precision and a tolerance of € = 10-2! the solution is converged to 
4 decimal places, which would have sufficed for this study. Regular 
double precision, however, would have been insufficiently precise 
and would not have led to the required accuracy. 

Although the orbits in figure 1 appear to follow a Keplerian orbit 
nicely, in figure 2 we show that this is not the case. Here we show 
a detail of the orbits of several stars from figure 1. Each orbital 
revolution is indicated with a thin line. The lines do not fully overlap, 
indicating that their orbits change with time. 

After having established a converged simulation for the S-stars (see 
figure 1), we introduce a relative shift by translating the Cartesian x 
coordinate of Star S5 by 107!0. This amounts to moving the initial 
position of S5 by ~ 15 m in the x direction. 

Instead of the separation in position space for each individual S- 
star, as illustrated in figure 1, we can also present the more usual 
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Figure 2. Magnification of a small section of the orbit of four S-stars over 
10 000yr of evolution. The orbit-to-orbit variation is small, and not noticeable 
in figure 1, but in this magnification the effect of chaos is visible in the non- 
overlapping orbits. 


total phase-space distance as a function of time, which we define as 
(Dejonghe & Hut 1986) 


n 
In(d) = sin ` li -Frol + Pii- Fiol]. (6) 
i=0 
Here 0 and 1 refer to the original and shifted solution, respectively. In 
figure 3 we present the evolution of the phase-space distance between 
the two solutions. The overplotted dotted line (in red) corresponds to 
a Lyapunov time scale of tq ~ 450 yr. 
In figure 4 we present the evolution of the root-mean-square differ. 
ence in the semi-major axes between the two solutions for the S-stars. 
Here we define 


1 n 
õa = 7 2 ldi = aial. ©) 
i=0 


The events are visible at discrete timesyfollowed by a constant 
difference ôa. In figure 3 a particularly strong‘event is visible in the 
phase-space distance evolution around t = 2876 years. According 
to punctuated chaos the orbital _separation/changes abruptly during 
events. We identify several with, thescolors and the vertical dashed 
lines (see figure 4). We automatically recognize these events in the 
simulation data using thesPelt,algorithm (Killick et al. 2012), using 
the piecewise constant model. Both free parameters in the Pelt al- 
gorithms (the controls the minimum distance between change points 
and the grid of possiblé-change points) were set to unity. 

To further investigate which stars are involved, we apply the Pelt 
algorithm to individual stars. Three examples are presented in fig- 
ure 5..Hete thé”colors are meant to guide the eye and indicate the 
events already identified on the full set of 27 stars in figure 4. The 
vertical red-dashed lines are the events as identified using the Pelt 
algorithm for each individual star. 

In table 1 we present the times at which events were detected, 


102 4 
10-2 
10-4 
6 
107 
10710 
10713 } —— low-bandpass filtered 
ta = 450 yr 
0 2000 4000 6000 8000 10000 


t [yr] 


Figure 3. Evolution of the phase-space distance (eq.6) for the\S-stars as a 
function of time (black). Here we adopted the phase-space differenc@between 
two solutions (the canonical initial conditions and thednitial conditions with 
a displacement of star S5 by 10~!9 in x); both are calculated using € = 10774 
(Lw = 128 bits). Overplotted (in orange) is a filtered version of the data using 
a Butterworth (1930) low-bandpass filter of-order_1. The dotted line gives a 
least-squares fit to the simulation data and indicates a Lyapunov time scale of 
ta = 450 year (red, as indicated). 
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Figure 4. Growth of the root mean square difference in the semi-major axes, 
Oa (eq.7), between the two solutions. The events are identified with the 
vertical red dashes. The color changes are introduced when the events lead 
to a growth in a. Two negative events are also identified, one around 667 yr 
and at 3914 yr. The largest event occurs around 2876 years. 


which stars are associated with these events, and at what distance 
from the black hole (in terms of au and in distance with respect to 
the sum of the stars’ Hill radii). The stars S2 and S18, are frequent 
participants, and they form the main drivers of chaos. We also present 
in table 1 the values for r;;, and r;;/ry. The former is defined as the 
distance between the two stars i and j nearest to the black hole at 
the moment of closest approach to the black hole. The latter gives 
the relative distance in terms of the radius of the Hill sphere of the 
two participating stars (second column). The last column gives the 
relative velocity between the two S-stars at the moment of closest 
mutual approach. 

In figure 6, we present the results of applying the theory of punc- 
tuated chaos to the evolution of the S-star cluster over 10 000 years. 
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Table 1. Moments of close encounters between two of the S stars (T, second 
column) and the black hole. The third column gives the distance (D, rounded 
to au) from the closest of the two stars to the black hole. The subsequent 
column gives the mutual distance between the two S-stars (r; j), followed by 
the mutual distance (r;;) as fraction of the sum of the S-star’s Hill radii. The 
last column gives the relative velocity (rounded to 10km/s) between the two 
encountering stars (in km/s). Each of these events can also be identified in the 
growth of the semi-major axis of the two identified stars (see also figure 4). 
This data can be used directly in the small routine to calculate the phase-space 
distance evolution of the S-stars using our presented theory on Punctuated 
Chaos (see appendix). Note that in the big event at T = 2876 yr, also stars S9 
and S14 participate (see figure 7). 


T participating stars D rij vig /TH Vij 

[yr] [au] [au] [km/s] 
12.5 838,824 842 48.7 5.2 5060 
407.8 $24,838 952 43.3 4.6 4650 
807.3. S29,S21 373 49.5 6.5 6030 
1571.8  838,S2 287 32.6 8.1 5580 
1649.7  S2,S18 516 42.2 6.7 3350 
1665.7  S2,S38 282 38.0 7.0 5580 
2023.8 S5,S14 780 55.8 11.9 4800 
2745.9  S2,S31 448 52.1 15.8 6170 
2875.7  S6,S21 635 13.5 1.5 3900 
3972.0  S5,S9 488 56.5 6.5 3660 
4875.7 S8,S27 720 54.0 6.0 3520 
6878.3 S$2,S18 512 30.4 4.0 3210 
7007.0 S2,S18 479 50.5 6.7 5420 


Here we overplotted the actual data and the filtered curve, as in fig- 
ure 3. The red curve results from measuring the time and magnitude 
of the events in semi-major axis (see figure 4 and table 1), which 
are then applied directly to our model for punctuated chaos (see 
section 2.2). 

In the appendix we present the algorithm, based on equation 6 that 
is used to draw the red curve in figure 6. The events (moments of 
close encounter and the mutual distance between the two S-stars at 
peribothron) are measured in the N-body simulations, and presented 
in table 1. Note that several punctuated events listed in table 1 are not 
directly related to the jumps in semi-major axis presented in figure 4 
(see also equation 7). For individual stars these jumps are visible) 
but they tend to be obscured by the ensemble averaging. Each’of the 
events listed in table 1, however, do relate to a jump in semi+major 
axis in the S-stars associated with the close encounter,Ahese jumps 
are visible in the evolution of the semi-major axis for individual S- 
stars; as examples we present this evolution for S6,.S27sand S66 in 
figure 5, which also show the early (13 yr, 408-yr and 807 yr) and late 
(7007 yr) jumps listed in table | (red dashed lines), 

The comparison between the theory and the’simulations is striking. 
Although the red curve in figure 6 does not perfectly match the filtered 
curve (black), it is astonishingly similax, A slight deviation where the 
red curve overshoots the black ¢wrye, around t = 7007 yr, seem to 
be caused by several closeencounters between S2 and S18. In our 
model, we anticipate on _ stan alðĥñe encounters, and the effect of 
multiple encounters in short succession between the same stars can 
also cause the evolutionof to decay. 


3.2 The big event at 2876 years 


The largest event happens around £ = 2876 yr after the start of the 
simulations: 

In figure 7 we present the distance of the nearest star to the black 
hole (blue), and the distance between this star and its next nearest 
stellar neighbor (orange), around the moment of the big event at 
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t = 2876 yr. Along the top we identify the stars involved in their 
close proximity. It starts at t = 2870 with S2 being the closest to the 
black hole, and $14 being closest to S2. At t ~ 2872 star S6 takes 
over the closest position to S2 from S14. Att ~ 2873 star S6 becomes 
the closest star to the black hole and at t ~ 2874 star S21 becomes 
the closest neighbor of star S6, until a little before t = 2876 star S6 
and S21 have their mutual close encounter with each other and with 
the black hole. 

In figure 8 we present for the stars S6 and S21 the evolution of 
ôr (left panels), i.e. the spatial separation between the two orbits; 
the semi-major axis (middle) and the difference in orbital frequency 
(right panels) around the moment of the big event at t = 2876 yr. 
In the evolution of the spatial separation (left) the consequences of 
the event are visible. The main driver of this evolution is the almost 
instantaneous change in orbital energy (middle panel, expressed here 
in semi-major axis) and the separation in orbital frequency“(right- 
most panel). 


4 DYNAMICS OF PUNCTUATED CHAOSAN THE'S-STARS 


The picture of punctuated chaos in Sec.2.2 lackspany dynamical 
characterisation of the events. They are-simply circumstances which 
change the frequency of the orbital motion. Even at the close of 
Sect.3.1, in the construction of figure.6, the data on the changes in 
frequency is simply read fron’ ayfull numerical simulation. In the 
present subsection we attempt to,explain, in order of magnitude, how 
these changes arise, on the basis"f few-body interactions. 

For the purposes ofthis discussion, let us reduce the problem to 
a 3-body system consistingyof a black hole (mass M, particle 0) and 
two S-stars (mass w, particles 1,2). Let r; be the position vector of 
star i relativesto the black hole, and focus on star 1. Its equation of 
motion is 

G(M+m)r,; Gm(ri-r2) Gmr 
~ 3 E EE =B (8) 
1 12 2 


fy 


r r; 


with obvious notation r;, r12. The first term on the right is the accel- 
eration of the motion of star 1 relative to the black hole, the second 
térm is the direct perturbation by particle 2, and the third term is the 
indirect perturbation, caused by the acceleration of the black hole by 
particle 2. (If we were to write down the equation for the case of N 
S-stars, the second and third terms on the right would be replaced by 
sums from j = 2 to N, with the index 2 replaced by j.) 


4.1 Direct perturbations 


From the previous section (especially figure 7 and table 1), we know 
empirically that the major events are associated with a close approach 
(to a distance of order 10—50 au) between two stars (see table 1), at 
times where the distances to the black hole are larger by an order of 
magnitude. If we label the stars in the encounter as 1 and 2 (ignoring 
all other stars for the time being), then it is clear that the perturbation 
on the motion of particle 1 is dominated by the second term on the 
right of eq.(8), i.e. the direct perturbation, due to the gravitational 
attraction of star 2. For the moment we also ignore the third term, 
but return to it briefly in Sec.4.2. 

We now consider the growth of chaos in the single S-star which we 
have hitherto labelled as star 1, and suppose that it has a succession 
of close encounters with other stars (in “events”) at some interval 
At. Let 6x denote the variation in some quantity x, such as energy or 
position, between two neighbouring solutions. Specifically, let 6Eg 
be the variation in specific energy E = |v?/2 —GM/r| at a time just 
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before an event, where v is the speed, and let ôrọ be the variation in 
position at the same time. (Henceforth the subscripts denote an index 
of events in a sequence of events.) 

Then we have 


6E; = ôEọ + 6AE, (9) 


where SE is the variation just before the next event, which will be 
the same as the value just after the first event; and AE is the change 
in E at that first event. Now we estimate this as 


AE = väv ~ VGM fr S, (10) 


where d is the distance of closest approach in the encounter, and V 
is the relative speed of the two stars in the encounter. Both these 
parameters are listed for the punctuations measured in the S-star 
cluster in table 1. (This estimate is made by taking Av to be the 
maximum perturbation Gm/d? multiplied by the time scale of the 
encounter, d/V.) We use the same estimate yGM /r for V as for v, 
so that 


AE ~ Gm/d, (11) 
Its variation can be estimated as 
6(AE) ~ St oro, (12) 
and so eq.(9) becomes 
OE, ~ Eo + OM Sri: (13) 
d2 


After this event the variation in the orbital frequency will be 6w , 
and so 


or, ~ ôro + wy aAt, (14) 


where a is the semi-major axis, just as in the simplified model of 
punctuated chaos described in Sec.2.2. Since wa? = GM and E ~ 
GM /a (though strictly here, M should be replaced by M +m, but we 
ignore the stellar mass compared with the black hole mass.), we can 
reexpress this as 


ôrı ~ ôro + OE; At/VGM/a. (15) 


By eq.(13), this in turn becomes 


ôrı ~ ôro + (6Eo + S oro) VGN] 016) 


Equations (13) and (16) are explicit estimatesewhich‘allow us to 
map the effect of an encounter on E and ôr. It.hasvmatrix 


A= 1 Gimhjd2 
~\At/ VGM/a  1+(Gm/d2)At/GM/a]} ` 


The larger eigenvalue of this matrix gives’ the factor by which the 
variation, in either E or r, is multiplied ås a result of an event. 

To get a feel for the magnhitude‘of this evolution, consider the six 
encounters in Table 1 which involve star S2. The mean distance of 
closest approach in these ‘is about 40 au, which we take for the value 
of d. Since the durationof the simulation is 10000 yr, for the mean 
time between these six events we adopt At = 1700 yr. The semi- 
major axis Of,S2us_as2 ~ 1000 au, and we take m = 20Mo and 
M = 4x‘T0° Me,as we have done throughout. (These units imply 
that G 4n*~r40.) Then the quantity 


Gm/d°At/VGM/a ~ 2, (17) 


and the larger eigenvalue 4 = 4. Thus the variations increase by this 
factor in 1700 yr, or by almost 4 dex in 10+ yr. While this estimate is 


only about half of the logarithmic growth seen in figures 3 and 6, it 
only includes half of the events listed in Table 1. The remaining events 
involve stars other than S2, which is why they were omitted from this 
estimate, but we show in Sec. 4.2 how their influence spreads to all 
stars, S2 included. 

The foregoing argument can be criticised on several grounds, and 
the main one is that it still depends on the numerical integration, as it 
draws on the distance of the observed two-body encounters. This may 
be estimated independently as follows. From the initial conditions of 
the integration we can estimate the central number-density n of the S- 
star cluster from the distance of the 6th-nearest neighbour to the black 
hole (Hénon 1971), which is about 2000 au. Hence n = 2 x 107!°, 
The median semi-major axis a is about 3000 au, and so we may 
estimate the typical speed v = ~GM/a as about 200 au/yr. It follows 
that a typical S-star 7 will experience one two-body encounter, with 
another S-star j in 10+ yr at a distance less than about rj; + 60 au. 
One or two will do so at still smaller distances, consistentavith what 
is seen in the numerical results. 


4.2 The role of the indirect acceleration 


While Sec.4.1 makes a reasonable case for supposing that direct 
encounters between pairs of S-stars aré*‘a major/driver of chaos in 
the system, it is not clear how many*stars)are affected, as only a 
fraction of all S-stars appear in Tables‘ Of course all stars are 
subject to two-body direct perturbations, but it might be thought that 
weak perturbations would drive chaos on a longer timescale. We now 
consider the effect of indirect perturbations, and will show that these 
affect all stars. Furthetmore they keep the growth of chaos in all stars 
in lockstep. 

The indirect perturbation of star 2 in eq.(8) is given by the third 
term on the tight,\and clearly depends on its distance r2 from the 
black hole. Thus we shall consider the pericentre passages of star 
2 asthe driver of this component of punctuated chaos. In one such 
eyent,\whose duration is of order q/v2, where q is the pericentre 
distance of star 2, the change in energy of star 1 may be estimated as 


AE ~ vjGm/(qv2). (18) 


Much as in Sec.4.1, we estimate vy ~ ~GM/q, but for the first star 


we take vı ~ GM /a, where a is the semi-major axis of star 1, since 
this star can be anywhere on its orbit when the event occurs. Thus 


AE ~ Gm/¥Vaq, (19) 


and so the variation in this quantity is 


(AE) ~ 6qgGm/Vaq?. (20) 


Now suppose that star 2 is one of those which is vigorously af- 
fected by direct interactions with other stars, and its variations have 
Lyapunov exponent 4. Then we can estimate that 


ôq ~ qo exp(at), (21) 


where ôqọo is some constant. Next, substituting this estimate into 
eq.(20), and adding the effect of a succession of such events, we 
readily see that the sum of the corresponding exponentials can be 
estimated by the effect of the most recent event. Thus the variation 
in the energy of star 1 at time t may be estimated as 


SE ~ ôqo exp(At)Gm/,/aq?. (22) 


Thus we conclude that chaos in star 1 grows with the same Lyapunov 
exponent as star 2, and this would be true even if star 1 were not 
subject to any direct perturbations. 
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In figure 9, we plot the time evolution of the separation in position 
space for each of the 27 S-stars, as well as for the central black hole. 
We observe that, on average, all particles diverge exponentially and 
with roughly the same rate. By fitting linear slopes to each curve in 
the figure, we measure an average individual Lyapunov time scale 
of 462 + 74 years. Figure 9 shows that all stars, and also the central 
black hole, diverge at the same rate on average. For the black hole 
the amplitude is smaller than for the bulk of the S-stars by roughly 
the ratio of their masses, i.e. 2 x 10°. 


5 CONCLUSIONS 


In order to explain short time-scale chaos in self-gravitating demo- 
cratic dynamical systems, we outline the theory of punctuated chaos. 
The essence of punctuated chaos is that instantaneous close encoun- 
ters, or events, drive the exponential growth due to ephemeral per- 
turbations, each of which is followed by a gradual drift in phase 
space. Our description of punctuated chaos delivers a qualitative and 
quantitative description of chaos in self-gravitating systems. 

We test this description on converged direct N-body simulations 
of self-gravitating systems under Newton’s equations of motion. We 
confirm the working of punctuated chaos on the chaotic orbits of the 
S-stars in orbit around the supermassive black hole in the Galactic 
center. The rather similar values of the mean orbital period of the 
S-stars and their Lyapunov time scale then is no coincidence. 

It turns out that also comet Halley has a Lyapunov time scale quite 
comparable to its orbital period (Boekholt et al. 2016). In this case 
the chaos is driven by interaction with Venus and Jupiter. We argue 
that the Lyapunov time scale under punctuated chaos is proportional 
to the mean interacting orbital period of the system, which for our 
selected systems is of the order of a few 100 years. 

According to our arguments and findings (Section 4.2), major 
bodies are affected on a similar time scale, but the amplitude 6 scales 
roughly with the ratio of the perturber to the perturbed mass, as one 
sees for the black hole in Fig.9. 

For the 27 S-stars in the Galactic center, we derive a Lyapunov 
time scale of 462 + 74 yr. This is somewhat comparable to the mean 
orbital period of the S-stars, or (Porh) = 269 +383 yr, consistent with 
the prediction from punctuated chaos. 

We also qualitatively compare the measured phase-space distance 
evolution with the theory of punctuated chaos (see Section,3° 1, es- 
pecially Fig. 6). 

In the theory (Sec. 3), the number of interactions\events) and the 
strength of these interactions are free parameters. In the comparison 
presented in Fig. 6, the moment and strength ofeach interaction are 
measured through an N-body simulation. But we also show (Sec. 4.1) 
that these parameters can be estimated imorđer of magnitude without 
resort to an N-body simulation, but from estimates of the density and 
kinematics of the inner S-stars. This procedure allows us to directly 
compare the theory withethe,simulations, as well as to identify the 
events that drive the exponential growth. 

If the previous/nalysis,holds, chaos manifests itself from repeated 
small perturbations that each induce a linear response in the separa- 
tion between neighbouring solutions of the system. The composition 
of these” linear fesponses resembles an exponential behavior and 
drives.chaos ifthe system. The perturbation, initiated by a close en- 
counter between two (or more) stars, is subsequently communicated 
to the rest of the system by the perturbed phase-space characteristics 
of the massive central body. In the case of the S-stars this is the 
super-massive black hole, and it is the Sun for comet Halley and the 
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rest of the Solar system. In this way the entire system is subject to 
the same smallest Lyapunov time scale. 
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APPENDIX: CALCULATING THE EVOLUTION OF THE 
PHASE-SPACE DISTANCE 


We calculate the evolution of the phase space distance, as presented 
in figure 6 (red curve) using the algorithms presented in the form 
of a Python script using the units module from the Astrophysics 
Multipurpose software Environment (AMUSE Portegies Zwart & 
McMillan 2018) in the listing below. The events (time 7 and mutual 
distance r;; between the two closest S-stars i and j near the black 
hole) used to draw the curves are presented in table 1. These events 
were automatically detected by identifying the closest encounters 
between at least two stars near the black hole. Here we included only 
those encounters with a mutual distance rj; < 60au. We introduce 
the constant k = 0.5 au which relates the effect of the perturbation to 
the orbital deviation. Looping over time, we calculate the new phase- 
space distance between the two solutions, given the perturbation 
introduced in the last close encounter (see equation 3). We empirially 
determine w œ (k/r; D from the N-body simulations by measuring 
rij and fitting k based on the evolution of the growth of the phase— 
space distance 6 (see equation. 6). The stated dependence on r;; here 
is justified by the d-dependence in equation 13. 


from amuse.lab import units 

def phase_space_distance(T, rij, 
delta = [0] | units.m 
t = [0] | units.yr 
dd = [1] | units.m 


k=4|units.au): 
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while t[—1] < T[—1]: 

d = dd[—1]*(1+((t{[—1]—T[len (dd) —1]). 
value_in(units.yr))*(k/rij [len(dd) —1]) 
+2) 

if t[—1] >= T[len(dd) ]: 

dd. append(d) 
t.append(t[—1l]+dt) 
delta . append (d) 

return t, delta 


The input arrays, T and rij are determined from close encounters 
between the S-stars and the black hole, as explained in Section 3.1, 
and present in table 1. Note that these encounters are not all the 
same as those identified using the Pelt algorithm in figure 4 because 
some close encounters do not show up as prominently in the global 
evolution of the orbital separation (see equation 7), but they do appear 
in the individual evolution of the semi-major axis of the S-stars. For 
S6, S21 and S66 these data are presented in figure 5. In particular 
the earlier encounters at 13 yr and 408 yr, and the late encounter at 
9046 yr do not appear from the analysis of the global change in semi- 
major axis through the Pelt algorithm, but they are visible as jumps 
in the evolution of the individual S-stars. 
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Figure 5. Evolution of the difference in semi-major axis’between the two con- 
verged solutions, for individual stars. From top to bottom, we show the results 
for stars S6, S27 and S66, respectively. The,vertical color bars are identical 
to those in figure 4. The vertical red dashed lines indicate the locations where 
the Pelt algorithm identifies theChange‘in a as events for the individual stars. 
The big event at t = 2876 years is visible in S6 and S66, but for S27 the event 
is less pronounced. A slightly earlier event, at 2771 years (this event is not 
listed in table 1, but prominentiin the right-hand panel in figure 8 for S21), is 
quite pronounced in $27, though. (It is to the left of the blue bar, whereas for 
S6 the event happens to the right side of the blue bar). 
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Figure 6. Normalized evolution of the phase-space distance, (eq.6) for the 
S-stars as a function of time (gray) with in black the¢filtered version of the 
data (see also figure 3). Overplotted is the result of.punctuatéd’chaos with 13 
close encounters (listed in Table 1) within rj; Æ 60 au between two of the 
S-stars i and j at a distance D < 1000 au fromthe, black hole (red). The green 
curve shows the result of ignoring all events, when the separation of the two 
orbits is driven by the initial displacement of stay S5 by 15 m in the positive 
x-direction (cf.eq.2). 
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Figure 7. Closest distance between stars as a function of time from t = 2870 yr 
to t = 2880 yr. The blue curve gives the distance from the black hole to the 
nearest star (identified in blue at the top). The orange curve gives the separation 
between this nearest neighbor and the next-nearest star (identified at the top 
row on orange). The vertical dotted lines indicate when the nearest neighbor 
changes (blue) or the next nearest neighbor (orange). The closest distance 
between S6 and S21 is reached at ~ 2876 yr at a mutual distance of ~ 13.5 au 
within 635 au of the supermassive black hole. 
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Figure 8. Illustration of an event driving the growth in separation between two solutions. We focus on the largest eventiat 2876\years, which was triggered by a 
close encounter between stars S6 (top row) and S21 (bottom row). From the simulation data, we plot the separation’ in. position space (left column), semi-major 
axis of one of the two solutions (middle), and the absolute value of the difference in orbital frequency betweensthem (right). We estimate mean values of a 
and w both before and after the event, which are sharply divided by the big event. Using these mean values, and the punctuated model for divergence from 
Sec. 2.2, we derive the analytical models for the evolution of ô, (left column). Furthermore, the energy exchange between S6 and S21 is conservative (to within 


a percent). Notice that the bottom right figre for S21 shows both events: at 2771 years and 2876 years¢ 
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Figure 9. The time evolution of the separation in position space (6,-) for 
each individual body, including the central black hole. On average, each 
body diverges at the same,exponential rate. Note, however, that these curves 
were calculated byscomparing the data from the two runs with € = 1077! 
and € =AQ-7*\ather than the perturbed and unperturbed solutions for 
e = 107%. 
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